home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / public / fax / src / sgi2fax / izoom.c < prev    next >
C/C++ Source or Header  |  1994-08-01  |  14KB  |  568 lines
  1. /*    $Header: /usr/people/sam/fax/sgi2fax/RCS/izoom.c,v 1.7 1994/02/28 14:23:10 sam Rel $
  2. /*
  3.  * Copyright (c) 1990, 1991, 1992, 1993, 1994 Sam Leffler
  4.  * Copyright (c) 1991, 1992, 1993, 1994 Silicon Graphics, Inc.
  5.  *
  6.  * Permission to use, copy, modify, distribute, and sell this software and 
  7.  * its documentation for any purpose is hereby granted without fee, provided
  8.  * that (i) the above copyright notices and this permission notice appear in
  9.  * all copies of the software and related documentation, and (ii) the names of
  10.  * Sam Leffler and Silicon Graphics may not be used in any advertising or
  11.  * publicity relating to the software without the specific, prior written
  12.  * permission of Sam Leffler and Silicon Graphics.
  13.  * 
  14.  * THE SOFTWARE IS PROVIDED "AS-IS" AND WITHOUT WARRANTY OF ANY KIND, 
  15.  * EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY 
  16.  * WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.  
  17.  * 
  18.  * IN NO EVENT SHALL SAM LEFFLER OR SILICON GRAPHICS BE LIABLE FOR
  19.  * ANY SPECIAL, INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY KIND,
  20.  * OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
  21.  * WHETHER OR NOT ADVISED OF THE POSSIBILITY OF DAMAGE, AND ON ANY THEORY OF 
  22.  * LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE 
  23.  * OF THIS SOFTWARE.
  24.  */
  25. /*
  26.  *      izoom- 
  27.  *              Magnify or minify a picture with or without filtering.  The 
  28.  *      filtered method is one pass, uses 2-d convolution, and is optimized 
  29.  *      by integer arithmetic and precomputation of filter coeffs.
  30.  *
  31.  *                              Paul Haeberli - 1988
  32.  */
  33. #include "stdio.h"
  34. #include "math.h"
  35. #include "assert.h"
  36. #include "izoom.h"
  37.  
  38. #define SHIFT           12
  39. #define ONE             (1<<SHIFT)
  40. #define EPSILON         0.0001
  41.  
  42. static FILTER *makefilt();
  43. static float filterrad();
  44. static float filterinteg();
  45. static float mitchell();
  46. static int (*xfiltfunc)(); 
  47.  
  48. static makexmap();
  49. static xscalebuf();
  50. static addrow();
  51. static divrow();
  52. static freefilt();
  53. static applyxfilt();
  54.  
  55. static int filtershape;
  56. static float blurfactor;
  57. static mults;
  58.  
  59. #define GRIDTOFLOAT(pos,n)      (((pos)+0.5)/(n))
  60. #define FLOATTOGRID(pos,n)      ((pos)*(n))
  61.  
  62. copyimage(getfunc,putfunc,nx,ny)
  63. int (*getfunc)(), (*putfunc)();
  64. int nx ,ny;
  65. {
  66.     int y;
  67.     short *abuf;
  68.  
  69.     abuf = (short *)malloc(nx*sizeof(short));
  70.     for(y=0; y<ny; y++) {
  71.         getfunc(abuf,y);
  72.         putfunc(abuf,y);
  73.     }
  74.     free(abuf);
  75. }
  76.  
  77. /*
  78.  *      general zoom follows
  79.  *
  80.  */
  81. zoom *newzoom(getfunc,anx,any,bnx,bny,filttype,blur)
  82. int (*getfunc)();
  83. int anx,any,bnx,bny;
  84. int filttype;
  85. float blur;
  86. {
  87.     zoom *z;
  88.     int i;
  89.     int xmults, ymults;
  90.  
  91.     z = (zoom *)malloc(sizeof(zoom));
  92.     z->getfunc = getfunc;
  93.     z->abuf = (short *)malloc(anx*sizeof(short));
  94.     z->bbuf = (short *)malloc(bnx*sizeof(short));
  95.     z->anx = anx;
  96.     z->any = any;
  97.     z->bnx = bnx;
  98.     z->bny = bny;
  99.     z->curay = -1;
  100.     z->y = 0;
  101.     z->type = filttype;
  102.     if(filttype == IMPULSE) {
  103.         if(z->anx != z->bnx) {
  104.             z->xmap = (short **)malloc(z->bnx*sizeof(short *));
  105.             makexmap(z->abuf,z->xmap,z->anx,z->bnx);
  106.         }
  107.     } else {
  108.         filtershape = filttype;
  109.         blurfactor = blur;
  110.         if(filtershape == MITCHELL) 
  111.             z->clamp = 1;
  112.         else
  113.             z->clamp = 0;
  114.         z->tbuf = (short *)malloc(bnx*sizeof(short));
  115.         z->xfilt = makefilt(z->abuf,anx,bnx,&z->nrows);
  116.         xmults = any*mults;
  117.         z->yfilt = makefilt(0,any,bny,&z->nrows);
  118.         ymults = bnx*mults;
  119. #ifdef notdef
  120.         fprintf(stderr,"TOTAL MULTS: %d\n",xmults+ymults);
  121. #endif
  122.         z->filtrows = (short **)malloc(z->nrows * sizeof(short *));
  123.         for(i=0; i<z->nrows; i++) 
  124.             z->filtrows[i] = (short *)malloc(z->bnx*sizeof(short));
  125.         z->accrow = (int *)malloc(z->bnx*sizeof(int));
  126.         z->ay = 0;
  127.     }
  128.     return z;
  129. }
  130.  
  131. getzoomrow(z,buf,y)
  132. zoom *z;
  133. short *buf;
  134. int y;
  135. {
  136.     float fy;
  137.     int ay;
  138.     FILTER *f;
  139.     int i, max;
  140.     short *row;
  141.  
  142.     if(y==0) {
  143.         z->curay = -1;
  144.         z->y = 0;
  145.         z->ay = 0;
  146.     }
  147.     if(z->type == IMPULSE) {
  148.         fy = GRIDTOFLOAT(z->y,z->bny);
  149.         ay = FLOATTOGRID(fy,z->any);
  150.         if(z->anx == z->bnx) {
  151.             if(z->curay != ay) {
  152.                 z->getfunc(z->abuf,ay);
  153.                 z->curay = ay;
  154.                 if(xfiltfunc) 
  155.                     xfiltfunc(z->abuf,z->bnx);
  156.             }
  157.             memcpy(buf,z->abuf,z->bnx*sizeof(short)); 
  158.         } else {
  159.             if(z->curay != ay) {
  160.                 z->getfunc(z->abuf,ay);
  161.                 xscalebuf(z->xmap,z->bbuf,z->bnx);
  162.                 z->curay = ay;
  163.                 if(xfiltfunc)
  164.                     xfiltfunc(z->bbuf,z->bnx);
  165.             }
  166.             memcpy(buf,z->bbuf,z->bnx*sizeof(short)); 
  167.         }
  168.     } else if(z->any == 1 && z->bny == 1) {
  169.             z->getfunc(z->abuf,z->ay++);
  170.             applyxfilt(z->filtrows[0],z->xfilt,z->bnx);
  171.             if(xfiltfunc)
  172.                 xfiltfunc(z->filtrows[0],z->bnx);
  173.             if(z->clamp) {
  174.                 clamprow(z->filtrows[0],z->tbuf,z->bnx);
  175.                 memcpy(buf,z->tbuf,z->bnx*sizeof(short)); 
  176.             } else {
  177.                 memcpy(buf,z->filtrows[0],z->bnx*sizeof(short)); 
  178.             }
  179.     } else {
  180.         f = z->yfilt+z->y;
  181.         max = ((int)f->dat)/sizeof(short)+(f->n-1);
  182.         while(z->ay<=max) {
  183.             z->getfunc(z->abuf,z->ay++);
  184.             row = z->filtrows[0];
  185.             for(i=0; i<(z->nrows-1); i++) 
  186.                 z->filtrows[i] = z->filtrows[i+1];
  187.             z->filtrows[z->nrows-1] = row;
  188.             applyxfilt(z->filtrows[z->nrows-1],z->xfilt,z->bnx);
  189.             if(xfiltfunc)
  190.                 xfiltfunc(z->filtrows[z->nrows-1],z->bnx);
  191.         }
  192.         if(f->n == 1) {
  193.             if(z->clamp) {
  194.                 clamprow(z->filtrows[z->nrows-1],z->tbuf,z->bnx);
  195.                 memcpy(buf,z->tbuf,z->bnx*sizeof(short)); 
  196.             } else {
  197.                 memcpy(buf,z->filtrows[z->nrows-1],z->bnx*sizeof(short)); 
  198.             }
  199.         } else {
  200.             memset(z->accrow,0,z->bnx*sizeof(int));
  201.             for(i=0; i<f->n; i++) 
  202.                 addrow(z->accrow, z->filtrows[i+(z->nrows-1)-(f->n-1)],
  203.                                                           f->w[i],z->bnx);
  204.             divrow(z->accrow,z->bbuf,f->totw,z->bnx);
  205.             if(z->clamp) {
  206.                 clamprow(z->bbuf,z->tbuf,z->bnx);
  207.                 memcpy(buf,z->tbuf,z->bnx*sizeof(short)); 
  208.             } else {
  209.                 memcpy(buf,z->bbuf,z->bnx*sizeof(short)); 
  210.             }
  211.         }
  212.     }
  213.     z->y++;
  214. }
  215.  
  216. freezoom(z)
  217. zoom *z;
  218. {
  219.     int i;
  220.  
  221.     if(z->type == IMPULSE) {
  222.         if(z->anx != z->bnx)
  223.             free(z->xmap);
  224.     } else {
  225.         freefilt(z->xfilt,z->bnx);
  226.         freefilt(z->yfilt,z->bny);
  227.         free(z->tbuf);
  228.         for(i=0; i<z->nrows; i++)
  229.             free(z->filtrows[i]);
  230.         free(z->filtrows);
  231.         free(z->accrow);
  232.     }
  233.     free(z->abuf);
  234.     free(z->bbuf);
  235.     free(z);
  236.  
  237. }
  238.  
  239. filterzoom(getfunc,putfunc,anx,any,bnx,bny,filttype,blur)
  240. int (*getfunc)(), (*putfunc)();
  241. int anx, any;
  242. int bnx, bny;
  243. int filttype;
  244. float blur;
  245. {
  246.     zoom *z;
  247.     int y;
  248.     short *buf;
  249.  
  250.     buf = (short *)malloc(bnx*sizeof(short));
  251.     z = newzoom(getfunc,anx,any,bnx,bny,filttype,blur);
  252.     for(y=0; y<bny; y++) {
  253.         getzoomrow(z,buf,y);
  254.         putfunc(buf,y);
  255.     }
  256.     freezoom(z);
  257.     free(buf);
  258. }
  259.  
  260. /*
  261.  *      impulse zoom utilities
  262.  *
  263.  */
  264. static makexmap(abuf,xmap,anx,bnx)
  265. short *abuf;
  266. short *xmap[];
  267. int anx, bnx;
  268. {
  269.     int x, ax;
  270.     float fx;
  271.  
  272.     for(x=0; x<bnx; x++) {
  273.        fx = GRIDTOFLOAT(x,bnx);
  274.        ax = FLOATTOGRID(fx,anx);
  275.        xmap[x] = abuf+ax;
  276.     }
  277. }
  278.  
  279. static xscalebuf(xmap,bbuf,bnx)
  280. short *xmap[];
  281. short *bbuf;
  282. int bnx;
  283. {
  284.     while(bnx>=8) {
  285.         bbuf[0] = *(xmap[0]);
  286.         bbuf[1] = *(xmap[1]);
  287.         bbuf[2] = *(xmap[2]);
  288.         bbuf[3] = *(xmap[3]);
  289.         bbuf[4] = *(xmap[4]);
  290.         bbuf[5] = *(xmap[5]);
  291.         bbuf[6] = *(xmap[6]);
  292.         bbuf[7] = *(xmap[7]);
  293.         bbuf += 8;
  294.         xmap += 8;
  295.         bnx -= 8;
  296.     }
  297.     while(bnx--) 
  298.         *bbuf++ = *(*xmap++);
  299. }
  300.  
  301. zoomxfilt(filtfunc)
  302. int (*filtfunc)(); 
  303. {
  304.     xfiltfunc = filtfunc;
  305. }
  306.  
  307. /*
  308.  *      filter zoom utilities
  309.  *
  310.  */
  311. static addrow(iptr,sptr,w,n)
  312. int *iptr;
  313. short *sptr;
  314. int w, n;
  315. {
  316.     while(n>=8) {
  317.         iptr[0] += (w*sptr[0]);
  318.         iptr[1] += (w*sptr[1]);
  319.         iptr[2] += (w*sptr[2]);
  320.         iptr[3] += (w*sptr[3]);
  321.         iptr[4] += (w*sptr[4]);
  322.         iptr[5] += (w*sptr[5]);
  323.         iptr[6] += (w*sptr[6]);
  324.         iptr[7] += (w*sptr[7]);
  325.         iptr += 8;
  326.         sptr += 8;
  327.         n -= 8;
  328.     }
  329.     while(n--) 
  330.         *iptr++ += (w * *sptr++);
  331. }
  332.  
  333. static divrow(iptr,sptr,tot,n)
  334. int *iptr;
  335. short *sptr;
  336. int tot, n;
  337. {
  338.     while(n>=8) {
  339.         sptr[0] = iptr[0]/tot;
  340.         sptr[1] = iptr[1]/tot;
  341.         sptr[2] = iptr[2]/tot;
  342.         sptr[3] = iptr[3]/tot;
  343.         sptr[4] = iptr[4]/tot;
  344.         sptr[5] = iptr[5]/tot;
  345.         sptr[6] = iptr[6]/tot;
  346.         sptr[7] = iptr[7]/tot;
  347.         sptr += 8;
  348.         iptr += 8;
  349.         n -= 8;
  350.     }
  351.     while(n--)
  352.         *sptr++ = (*iptr++)/tot;
  353. }
  354.  
  355. static FILTER *makefilt(abuf,anx,bnx,maxn)
  356. short *abuf;
  357. int anx, bnx;
  358. int *maxn;
  359. {
  360.     FILTER *f, *filter;
  361.     int x, min, max, n, pos;
  362.     float bmin, bmax, bcent, brad;
  363.     float fmin, fmax, acent, arad;
  364.     int amin, amax;
  365.     float cover;
  366.  
  367.     f = filter = (FILTER *)malloc(bnx*sizeof(FILTER));
  368.     *maxn = 0;
  369.     mults = 0;
  370.     for(x=0; x<bnx; x++) {
  371.         if(bnx<anx) {
  372.             brad = filterrad()/bnx;
  373.             bcent = ((float)x+0.5)/bnx;
  374.             amin = floor((bcent-brad)*anx+EPSILON);
  375.             amax = floor((bcent+brad)*anx-EPSILON);
  376.             if(amin<0)
  377.                amin = 0;
  378.             if(amax>=anx)
  379.                amax = anx-1;
  380.             f->n = 1+amax-amin;
  381.             mults += f->n;
  382.             f->dat = abuf+amin;
  383.             f->w = (short *)malloc(f->n*sizeof(short));
  384.             f->totw = 0;
  385.             for(n=0; n<f->n; n++) {
  386.                 bmin = bnx*((((float)amin+n)/anx)-bcent);
  387.                 bmax = bnx*((((float)amin+n+1)/anx)-bcent);
  388.                 cover = filterinteg(bmax,1)-filterinteg(bmin,0);  
  389.                 f->w[n] = (ONE*cover)+0.5;
  390.                 f->totw += f->w[n];
  391.             }
  392.         } else {
  393.             arad = filterrad()/anx;
  394.             bmin = ((float)x)/bnx;
  395.             bmax = ((float)x+1.0)/bnx;
  396.             amin = floor((bmin-arad)*anx+0.5+EPSILON);
  397.             amax = floor((bmax+arad)*anx-0.5-EPSILON);
  398.             if(amin<0)
  399.                amin = 0;
  400.             if(amax>=anx)
  401.                amax = anx-1;
  402.             f->n = 1+amax-amin;
  403.             mults += f->n;
  404.             f->dat = abuf+amin;
  405.             f->w = (short *)malloc(f->n*sizeof(short));
  406.             f->totw = 0;
  407.             for(n=0; n<f->n; n++) {
  408.                 acent = (amin+n+0.5)/anx;
  409.                 fmin = anx*(bmin-acent);
  410.                 fmax = anx*(bmax-acent);
  411.                 cover = filterinteg(fmax,1)-filterinteg(fmin,0);  
  412.                 f->w[n] = (ONE*cover)+0.5;
  413.                 f->totw += f->w[n];
  414.             }
  415.         }
  416.         if(f->n>*maxn)
  417.             *maxn = f->n;
  418.         f++;
  419.     }
  420.     return filter;
  421. }
  422.  
  423. static freefilt(filt,n)
  424. FILTER *filt;
  425. int n;
  426. {
  427.     FILTER *f;
  428.  
  429.     f = filt;
  430.     while(n--) {
  431.         free(f->w);
  432.         f++;
  433.     }
  434.     free(filt);
  435. }
  436.  
  437. static applyxfilt(bbuf,xfilt,bnx)
  438. short *bbuf;
  439. FILTER *xfilt;
  440. int bnx;
  441. {
  442.     short *w;
  443.     short *dptr;
  444.     int n, val;
  445.  
  446.     while(bnx--) {
  447.         if((n=xfilt->n) == 1) {
  448.             *bbuf++ = *xfilt->dat;
  449.         } else {
  450.             w = xfilt->w;
  451.             dptr = xfilt->dat;
  452.             val = 0;
  453.             n = xfilt->n;
  454.             while(n--) 
  455.                  val += *w++ * *dptr++;
  456.             *bbuf++ = val / xfilt->totw;
  457.         }
  458.         xfilt++;
  459.     }
  460. }
  461.  
  462. static float filterrad()
  463. {
  464.     switch(filtershape) {
  465.         case BOX:
  466.             return 0.5*blurfactor;
  467.             break;
  468.         case TRIANGLE:
  469.             return 1.0*blurfactor;
  470.             break;
  471.         case QUADRATIC:
  472.             return 1.0*blurfactor;
  473.             break;
  474.         case MITCHELL:
  475.             return 2.0*blurfactor;
  476.             break;
  477.     }
  478. }
  479.  
  480. static float quadinteg(x)
  481. float x;
  482. {
  483.    if(x<-1.0)
  484.        return 0.0;
  485.    if(x<-0.5)
  486.        return 2.0*((1.0/3.0)*(x*x*x+1.0) + x*x + x);
  487.    else
  488.        return -(2.0/3.0)*x*x*x + x + 0.5;
  489. }
  490.  
  491. static float filterinteg(x,side)
  492. float x;
  493. int side;
  494. {
  495.     float val;
  496.  
  497.     x = x/blurfactor;
  498.     switch(filtershape) {
  499.         case BOX:
  500.             if(x<-0.5)
  501.                 return 0.0;
  502.             else if(x>0.5)
  503.                 return 1.0;
  504.             else
  505.                 return x+0.5;
  506.             break;
  507.         case TRIANGLE:
  508.             if(x<-1.0)
  509.                 return 0.0;
  510.             else if(x>1.0)
  511.                 return 1.0;
  512.             else if(x<0.0) {
  513.                 val = x+1.0;
  514.                 return 0.5*val*val;
  515.             } else {
  516.                 val = 1.0-x;
  517.                 return 1.0-0.5*val*val;
  518.             }
  519.             break;
  520.         case QUADRATIC:
  521.             if(x<0.0)
  522.                 return quadinteg(x);
  523.             else
  524.                 return 1.0-quadinteg(-x);
  525.             break;
  526.         case MITCHELL:
  527.             if(side == 0)
  528.                 return 0.0;
  529.             return mitchell(x);
  530.     }
  531. }
  532.  
  533. static float p0, p2, p3, q0, q1, q2, q3;
  534.  
  535. /*
  536.  * see Mitchell&Netravali, "Reconstruction Filters in Computer Graphics",
  537.  * SIGGRAPH 88.  Mitchell code provided by Paul Heckbert.
  538.  */
  539. static mitchellinit(b,c)
  540. float b, c;
  541. {
  542.  
  543.     p0 = (  6. -  2.*b        ) / 6.;
  544.     p2 = (-18. + 12.*b +  6.*c) / 6.;
  545.     p3 = ( 12. -  9.*b -  6.*c) / 6.;
  546.     q0 = (           8.*b + 24.*c) / 6.;
  547.     q1 = (        - 12.*b - 48.*c) / 6.;
  548.     q2 = (           6.*b + 30.*c) / 6.;
  549.     q3 = (     -     b -  6.*c) / 6.;
  550. }
  551.  
  552. static float mitchell(x)        /* Mitchell & Netravali's two-param cubic */
  553. float x;
  554. {
  555.     static int firsted;
  556.  
  557.     if(!firsted) {
  558.         mitchellinit(1.0/3.0,1.0/3.0);
  559.         firsted = 1;
  560.     }
  561.     if (x<-2.) return 0.0;
  562.     if (x<-1.) return 2.0*(q0-x*(q1-x*(q2-x*q3)));
  563.     if (x<0.) return 2.0*(p0+x*x*(p2-x*p3));
  564.     if (x<1.) return 2.0*(p0+x*x*(p2+x*p3));
  565.     if (x<2.) return 2.0*(q0+x*(q1+x*(q2+x*q3)));
  566.     return 0.0;
  567. }
  568.